Spectral Analysis by the Method of Consistent Constraints 
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Two major challenges of numeric analytic continuation — restoring the spectral density, s(cj), from 
corresponding Matsubara correlator, g(r) — are (i) producing the most smooth/featureless answer 
for s(cj), without compromising the error bars on g(r), and (ii) quantifying possible deviations of 
the produced result from the actual answer. We introduce the method of consistent constraints that 
solves both problems. 



A dynamic linear-response function can be obtained 
from the associated Matsubara correlator by analytic 
continuation of the latter. A prototypical example is 
the ground-state single-particle Matsubara Green's func- 
tion in imaginary-time representation, g(r). If g(r) 
is specified numerically pQ, the procedure of analytic 
continuation — often referred to as spectral analysis — 
amounts to finding an appropriate spectral function s(u) 
related to g(r) by an integral equation. In our prototyp- 
ical case, the equation reads (here the function s(uj) is 
known to be identically zero at lj < and non-negative 
at u > 0) 



9(r) 



r s(uj) duo . 
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The notorious difficulty of the problem comes from the 
fact that finding s(u) does not reduce to the require- 
ment that the integral in the right-hand side of (fil) re- 
produces the values of g(r) within their error bars. Be- 
ing ill posed, i.e. subject to the sawtooth instability, the 
problem of numerically finding s(u) features infinitely 
many solutions, ranging continuously from very smooth 
to extremely noisy ones. It is crucial, thus, not only to 
satisfy Eq. ([!]) for a specified set of r-points, but also to 
guarantee that s(uj) is free of sawtooth artifacts. The 
following two approaches work rather well in achieving 
this goal: (i) the stochastic-optimization method (SOM) 
[2] and (ii) the maximum entropy method (MEM) [3j H] . 
Within SOM, one employs a stochastic process of mini- 
mizing the standard x 2 -measure to produces a large num- 
ber of noisy solutions, and then takes an average of all 
the solutions, which is a legitimate procedure thanks to 
linearity of Eq. (fTJ. In the statistical limit, the outcome 
of SOM is a rather smooth function as the sawtooth ar- 
tifacts are averaged out. By construction, the accuracy 
of reproducing g{r) with s(uo) is not compromised, but, 
speaking generally, with SOM one cannot guarantee that 
the final result for s(uo) is the smoothest of all the func- 
tions consistent within the error bars of g(r). MEM is 
complimentary to SOM in the sense that, on one hand, it 
does guarantee that the outcome for s(uj) is the best one 
within a certain class of smooth functions (this class is 
selected by formulating a "target" or "default" model for 
the solution, see [SJH]). On the other hand, smooth so- 
lutions are produced at the expense of a systematic bias 



introduced by the default model; the bias becomes more 
pronounced as the error bars on g(r) are decreased. 

Existing methods of spectral analysis, with MEM and 
SOM as characteristic examples, seem to follow the gen- 
eral principle (cf., e.g., Tikhonov-Phillips regularization 
methods [5-7 ) that there is always a compromise be- 
tween the requirement of s(uj) being as smooth as possi- 
ble and the requirement of reproducing g(r) within the 
error bars. 

In this Letter, we show that the "compromise prin- 
ciple" is a mere prejudice. It is possible, and relatively 
easy (!), to meet the condition of smoothest possible s(u) 
while perfectly respecting the error bars on g(r). The 
price that one has to pay for this luxury is the necessity to 
introduce a feedback loop locally adjusting the smooth- 
ness constraints on s(uj) to ensure consistency with the 
error bars on g(r). More importantly, the method of 
consistent constraints (MCC) has a simple built-in tool 
of quantifying the accuracy of s(u). 

The issue of quantifying the accuracy of s(lj) is yet an- 
other challenge for the problem of spectral analysis. The 
sawtooth instability implies that any error bar on s(uj) 
should necessarily be of a conditional character. The con- 
dition which we adopt (finding it natural) is as follows. 
Any function s(u) that we include into the class of le- 
gitimate deviations from the optimal (i.e. most smooth) 
solution is subject to the constraint of not having ex- 
tra qualitative features like, e.g., extra bumps or min- 
ima (more specifically, the constraint is to keep the same 
structure of sign-domains of the second derivative, see 
below). With this constraint — clearly justified in the 
asymptotic regime of appropriately small error bars on 
the function g(r) — we employ MCC to deliberately dis- 
tort the function s(uj) at a certain uj = uo* and see the 
extent to which the distortion remains consistent with 
the error bars of g(r). 

Illustrative examples. Leaving technical description of 
the MCC numeric protocol to the second half of the pa- 
per, here we consider two illustrative examples. Figures 
IT] and ^1 show the results (cases A and B, respectively) 
of applying MCC to determine spectral functions for 
g(r) specified numerically with known error bars. Corre- 
sponding functions s(uj) are shown along with their exact 
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FIG. 1. Case A: The function s(uo) (circles), obtained by 
spectral analysis, is very close to the exact spectral density 
s e (cj) (dashed line). The task for the error analysis is to 
confirm small error bars on s(uj) without knowing s e (uj). 



counterparts s e (uj). The functions 



g 3 (r) = / e UT s(uj) duj 
Jo 



(2) 



coincide with g(r) within the error bars on the latter, 
while the functions s(uj) are quite smooth. Now we have 
to characterize possible deviations of s(u) from s e (uu) — 
pretending that the latter is unknown — making sure that 
what we get is consistent with the deviations we see in 
the two figures. In Fig. [3] we show how we quantify the 
uncertainty of the sharpness of the second peak. The 
error bars for both cases are presented in Figs. [4] and [5) 
The asymmetry of the error bars reflects the tendency of 
the reference solution to broaden sharp features. Apart 
from the error bar asymmetry, correlations between the 
errors are very informative, as is clear from Fig. [3] 

Objective function. Numerically, the left hand side 
of Eq. ([I]) is known on a discrete set of r-points r = 
(ti,T2, . . . ,tjv) and gfc) values have statistical error 
bars Gi. The solution for the spectral function will 
be defined on a dense enough set of frequency points 
uj = (u;i,u;2, • • • ,^jv) so that the integral in Eq. is 
transformed into the finite-sum expression defining a set 
of g s (ji) values 
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(3) 
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The objective function to be minimized in the process 
of searching for the smooth solution s(ujj) involves several 
terms O = ^2 k 0&. In what follows we describe the min- 
imal number of objectives which allow one to achieve the 
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FIG. 2. Case B: Close to the second peak, the function 
s(uj) (circles), obtained by spectral analysis, is substantially 
smother than its exact counterpart s e (cj) (dashed line). The 
challenge for the error analysis here is to characterize possible 
deviations of s(oj) from s e (cj), without knowing the latter. 



final goal (all results presented in this Letter were based 
on them) . The first and most important term is the stan- 
dard x 2 measure which penalizes differences g(ri) —g s (ri) 
outside of the computed error bars, c^. For simplicity, we 
write this measure for the case of uncorrelated statistical 



o 1= iVx 2 = E 



r g(n)-9s(Tiy 2 



(4) 



The major goal is to have this objective of the order of 
unity. Penalty for first derivatives is preventing the devel- 
opment of the saw-instability, i.e. fast changing solutions 
are disfavored 



M-l 



(5) 



To simplify notations, we introduced dj = |«j+i — Sj\. 
Since the minimization of quadratic form does not guar- 
antee that the solution is non-negative, we also need a 
penalty which suppresses the amplitudes of the solution 
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(6) 



Finally, we introduce penalty for the solution to deviate 
from some "target" form Sj for j = (2, . . . , M — 1): 



M-l 
3=2 



(7) 




FIG. 3. Example of error analysis for the case B. Shown with 
a solid line is the reference (the smoothest) solution s(cj), pre- 
viously presented in Fig. [5] The solution s*(cj) is obtained by 
pulling the solution up — without compromising the deviation 
of g s (j) from s(ui) — at the point u — 2.5 corresponding to the 
maximum of the second peak. The solution s*(cj) corresponds 
to the threshold of appearance of an extraneous feature, a 
shoulder at uj w 6.5 (that will develop into a bump if we keep 
pulling up). In this sense, s*(cj) characterizes maximal po- 
tential deviation of the reference solution from the exact one 
in the vicinity of the second peak, cf. Fig. [2] 
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FIG. 4. Case A with three characteristic error bars (cf. 
Fig.Q. 
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There is nothing new in the idea of introducing regu- 
larization measures similar to O2 , or Q4 [3j EHZ] but in 
the past it was done with j-independent coefficients con- 
sidered essentially as input parameters (ultimately opti- 
mized to achieve the best, but somewhat compromised, 
X 2 ). In our scheme, it is absolutely crucial that con- 
straints control every point of the solution and are ad- 
justed by the feedback loop to be consistent with the 
properties of the solution itself. Only in this case do we 
have a guarantee that in the limit of vanishingly small 
error bars on g(r) the final solution will always reach the 
X 2 ~ 1 limit. 

MCC protocol Given an objective based on the 
positive definite quadratic form for Sj , one can relatively 

easily find the solution minimizing it, s op (O), e.g. by 
the gradient method. This solution, however, may have 
two serious drawbacks: It may contain negative values 
of Sj and be far from meeting the crucial requirement 
of having % ~ 1. The following self-consistent iterative 
protocol is designed to eliminate both shortcomings. 

1. Let index k denote the number of performed itera- 



tions. Start with k = 0, some initial solution Sj , large 

and zero values of 



penalties for first derivatives Dj , 
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FIG. 5. Case B with the error bars (cf. Figs. [5] and [3]). The 
error bars are essentially asymmetric, reflecting the tendency 
of the reference solution to broaden sharp features. 
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We find it sufficient to have T = D- for any k (with 
small modification in the protocol quantifying error 
bars, see below), but one is free to design other rules for 
these coefficients. 



2. Determine the next iteration solution by minimiz- 
ing the objective function, Sj = Sj^O^); set 
k->k + l. 



3. Since Q\ term is the only reason for having a non-flat 
result, the solution is now analyzed to adjust the objec- 
tive so that penalties compromising the x 2 measure are 
reduced, and penalties for developing the negative solu- 



.(fc-i) 



tion are increased: If d K - ' exceeds C /D K - L) (in practice, 
we use C — 0.1) then 



£><*> = C/df 



(8) 



Otherwise, the penalty is considered to be too conserva- 



i(fe) 



(k) 



tive and is increased as D K - = 2D y - . To prevent diver- 
gent behavior, determine £> m in = min{Dj} and restrict 
allowed values to r\D m i n , where r is some large number 
(typically of order of 10 3 ). After that set T J (/c) = Df\ 

(k) (k) 

4. If s- < introduce large penalty A y - J = A max (it 
can be as large as 10 8 ) to prevent the solution from 

(k) 

going negative in the next iteration. If s) > decrease 



the penalty, A 



(k) 



.(fc-i) 



/10 (for positive values of the 



solution A-coefflcients decay exponentially fast with the 
number of iterations) . 

dj — v 5 j+i + s j-i)/2 to stabilize second 
derivatives of the solution. This determines the new 



5. Finally, set sjW = (s!- 
luti 
objective function 0^ k \ Proceed to step 2. 



In essence, the procedure adjusts regularization param- 
eters by feedback from the solution itself and ultimately 
has them small enough to admit a solution of Eq.Q 
within the error bars, and large enough to have a smooth 
solution. If there are ^-functional peaks in the spectral 
function one should add them to the solution at some lo- 
cations and exclude their amplitudes from the objectives 
O2 and O4 for obvious reasons. The nature and stability 
of the MCC feedback loop is also compatible with inter- 
rupting it from time to time and suggesting modifications 
to the existing solution which minimize only the y 2 — the 
subsequent MCC protocol quickly erases unwanted rip- 
ples. We also find it useful to "level-off" large point-to- 
point fluctuations in the regularization parameters after 
several feedback loops. 

The protocol of examining error bars on Sj is iden- 
tical to the one described above except that for some 
frequency point uj m = uj* the target parameter s m re- 
mains fixed at some predetermined value s* and T m is 
made large enough to suppress significant deviations of 
the solution from s*. As we increase/decrease the value 



of s* we ultimately observe that either the value of \ 2 
is increased by a factor of two, or a new feature appears 
on the curve. This determines the error margins on the 
final solution. [Since in the protocol of determining error 
bars on Sj the values of s m and T m are excluded from the 
feedback loop the x 2 measure may increase as we "pull" 
the point due to hysteresis effect.] 

Discussion. As is seen from the above-described pro- 
tocol, the central idea of MCC is to iteratively adjust — 
tighten/relax — regularization constraints, by feedback 
from the solution s(uj) that minimizes the objective func- 
tion for the previous iteration. The choice of constraints 
is rather wide; they can deal with the values of the func- 
tion s(cj), as well as with the values of its first and 
higher derivatives; within one and the same iteration, 
or addressing previous iterations as well. The crucial 
common feature is the locality of constraints in cj-space, 
which readily allows one to judge — see steps 3 and 4 
of the protocol — whether constraints are (i) too restric- 
tive, (ii) too loose, or (hi) perfectly consistent. Perfor- 
mancewise, it is due to those simple local measures of 
consistency of each specific constraint that feedback it- 
erations rapidly fine-tune the objective function to be 
the most restrictive — while remaining free of systematic 
bias — within the error bars of the function g(r). 

In principle, the set of consistent constraints can be 
extended to include any knowledge about likely features 
of the function s(uj). For example, expecting a simple 
monotonic asymptotic behavior in the tail, one may opt 
to severely penalize wrong signs of (higher) derivatives, 
similarly to penalizing negative values of s(cj). 

The MCC protocol can be applied to any problem of 
restoring a function from a representative set of its in- 
tegrals with an arbitrary kernel. An immediate example 
is the numeric function g(r) itself, when it comes from a 
Monte Carlo simulation in terms of integrals of g(r) over 
a number of intervals ("bins"). An accurate value of g{r) 
for any r can then be extracted by the MCC — either di- 
rectly, or as a by-product of restoring s(u) from the bin 
integrals. 

In absorption spectroscopy, MCC can be used for unbi- 
ased restoring unknown density distribution from binned 
absorption images, or just to produce a smooth — still un- 
biased (!) — image from data integrated over bins. 
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